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ABSTRACT 

We present a reconstructions of galaxy-cluster-scale mass distributions from simulated 
gravitational lensing data sets including strong lensing, weak lensing shear, and mea¬ 
surements of quadratic image distortions - flexion. The lensing data is constructed 
to make a direct comparison between mass reconstructions with and without flex¬ 
ion. We show that in the absence of flexion measurements, signiflcant galaxy-group 
scale substructure can remain undetected in the reconstructed mass profiles, and that 
the resulting profiles underestimate the aperture mass in the substructure regions by 
^ 25 — 40%. When flexion is included, subhaloes down to a mass of ^ 3 x 10^^ Mq can 
be detected at an angular resolution smaller than 10". Aperture masses from profiles 
reconstructed with flexion match the input distribution values to within an error of 
^ 13%, including both statistical error and scatter. This demonstrates the important 
constraint that flexion measurements place on substructure in galaxy clusters and its 
utility for producing high-fidelity mass reconstructions. 

Key words: gravitational lensing: strong, gravitational lensing: weak, galaxies: clus¬ 
ters: general, dark matter, methods: data analysis 


1 INTRODUCTION 

The distribution of mass structure in the Universe provides a 
fundamental test for any plausible cosmological model. Em¬ 
pirical data from multiple experiments have robustly con¬ 
strained cosmology to be consistent with a flat ACDM model 


([Hinshaw et al.||2013 Kowalski et al.||2008 Planck Collab- 


[oration et al. 2013 Eisenstein et al. 2005). But to further 
understand the constituents of this cosmology, in particular 
the dark matter, we must investigate the properties of the 
structures which form from the primordial mass distribu¬ 
tion at many scales. To understand the interaction between 
dark matter and baryonic matter, or to understand the self¬ 
interactions between dark matter particles we must mea¬ 
sure the distribution of structures on relatively small scales, 
meaning measuring structures down to galaxy scales where 
the properties of both baryons and dark matter influence the 
evolution of structure. As interactions between dark matter 
and baryons are more important to the structure of haloes 
for more non-linear overdensities, measurements of substruc¬ 
tures inside of larger halos (e.g., substructures within galaxy 
clusters) provide an important window into the nature of 
dark matter. 

Eor a ACDM cosmology, the particular initial condi- 
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tions (e.g., the primordial matter power spectrum) as well 
as the particular properties of the dark matter - both its 
self-interactions and its interaction with baryons, if any - 
will affect the structures that are formed. This includes the 


mass function of virialized dark matter haloes (Tinker et al. 


2008) and the spatial and mass distributions of subhaloes 
"(Nagai & Kravtsov 2005| Kravtsov et al. 2004). The exis¬ 
tence of roughly self-similar haloes at all scales is a robust 
prediction of ACDM ( Moore et al.|1999 ), making the detec¬ 
tion of galaxy cluster substructures an important cosmolog¬ 
ical test. The fraction of a cluster’s mass which is associated 
with substructures as well as the other subhalo properties 
(e.g., density profiles, ellipticities, etc.) can be compared to 
subhaloes produced in cosmological A-body simulations to 
further constrain dark matter properties. Recent cosmologi¬ 
cal simulations of clusters with self-interacting dark matter 
(SIDM) have found that dark matter interactions affect both 
the shape of haloes, with more spherical profiles and flatter 
halo cores as a result of a higher interaction cross-section 
for galaxy-scale haloes, as well as reducing the number of 
subhaloes of a given mass ( Shaw et al.||2006 Rocha et al 


2013 


Peter et ^|2Q13 ). 


To measure the mass distribution of substructures 
within galaxy clusters we cannot use statistical measures of 
cluster properties, such as empirical scaling relations (see, 
e.g., Giodini et al. 2013 for a review of these relations). 
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These relationships are useful for large samples of clusters 
but the intrinsic scatter between individual clusters limits 
their predictive power. Neither can we use baryonic trac¬ 
ers of structure due to uncertainty in the interactions be¬ 
tween baryons and dark matter. However, gravitational leas¬ 
ing, and particularly combining strong lensing deflections 
and linear image distortions (weak lensing shear and con¬ 
vergence) has been successfully shown to produce detailed 
mass reconstructions of galaxy clusters and resolved compli¬ 
cated, asymmetric morphologies (Clowe et al.|2006 Natara- 
jan et al.||2Q07| |2009| ). 


For typical cluster imaging datasets, the correlation of 
spatially-separated galaxy shapes with uncertain redshifts 
produces an effective smoothing length of > 10 '' away from 
the central region of the cluster which is constrained by 
strong lensing. Simulations indicate that at least half of dark 
matter subhaloes will exist outside of the Einstein radius of 
a the main halo ( Nagai Kravtsov|20Q5 ). In this large por¬ 
tion of the cluster, at least 99% of the projected area within 
the cluster virial radius, substructures cannot be efficiently 
detected in galaxy clusters without additional information 
or the investment of prohibitively large amounts of telescope 
time. 

In order to overcome these limitations, to improve the 
resolution of reconstructed mass distributions, and thus to 
probe small-scale mass structure, we includie higher-order 
lensing effects. In regions where the reduced lensing shear 
g — 7/(1 — /^) is a significant fraction of unity, as is the 
case within a radius of ~60" from the cluster center for a 
typical galaxy cluster, nonlinear image distortions become 
significant and the linear distortion approximation becomes 
increasingly less valid as the shear increases (see [Bartel- 


|mann Schneider|| 200 f] for a review of linear weak lens¬ 
ing), which also degrades the reliability of galaxy ellipticity 


as an estimator for shear with standard approaches (Mel- 
[chior et "^^12010 ). Quadratic image distortions due to non- 


negligible third derivatives of the lensing potential, called 
flexion, can be observed and quantitatively measured in 
this regime. Flexion distortions cause intrinsically elliptical 
galaxy images to be bent by into “banana-shaped” arclets, 
augmenting the image distortions due to linear shear and 
magnification. The two lensing fields which cause flexion dis¬ 
tortions are sourced by the gradients of the convergence and 
shear fields. This has the practical effect of making flexion a 
sensitive probe of small-scale variations in the surface mass 
density. Limits on the measurable flexion also constrain the 
amount of possible substructure in a region where no flexion 
is significantly detected. 

In this paper we present a method for reconstruct¬ 
ing a lensing mass distribution using multiple image posi¬ 
tions (strong lensing), ellipticity measurements (weak lens¬ 
ing shear), and flexion measurements. This significantly im¬ 
proves the sensitivity of the reconstruction to small mass 
substructures; we are able to resolve structures to a signif¬ 
icantly finer spatial resolution and a lower mass threshold 
than with strong lensing and shear alone. We are able to 
directly detect structures with mass ~ 9 X IO^^Mq which 
are not detected without the flexion measurements, even in 
the case where they are located far from the strong lensing 
regime. 

Throughout we employ a standard flat cosmological 
model with Vtm — 0.3, Qa = 0.7, and Hq — 70 km/s/Mpc. 


We also employ standard complex notation for lensing fields, 
with V = + 182 being the derivative operator acting in 

the image plane. 


2 MASS RECONSTRUCTION METHOD 

Our mass reconstruction method is an expansion of the 


method described in Bradac et al. ( 

2009), which itself was an 

improvement of previous methods ( 

Bradac et al. 2005 2006). 

We review the basic concept of the method here to place our 


expansion in context. We include five complex lensing fields, 
all of which yield measurable lensing distortions which are 
inputs into our reconstruction. All five fields arise as linear 
combinations of the first through third derivatives of a sin¬ 
gle, real-valued lensing potential ^ sourced by the lens mass 
distribution. These lensing fields are the deflection a = V^, 
the convergence | VV*'?/, the shear 7 = | I-flexion 
(comatic flexion) = V/^, and 3-flexion (trefoil flexion) 
0 = V 7 . Here we have used complex notation for both the 
fields and the derivative operators, and the star operator de¬ 
notes the complex conjugate. The positions of the multiple 
images in a strong lensing system constrain the deflection, 
and the ellipticity of lensed images of background galaxies 


constrains the shear and the convergence (see Bartelmann 


Schneider|| 20 ^ for a review of strong lensing and shear 
constraints on lens mass distributions). 

Both I-flexion and 3-flexion are constrained by mea¬ 
surements of non-linear distortions due to gradients in the 
convergence and shear. Since the convergence is the sur¬ 
face mass density of the lens scaled to the critical lensing 
density, measurements of image distortions due to flexion 
provide an important direct probe of mass gradients, and 
therefore small-scale mass structures. The flexion fields can 
be measured using several methods available in the litera- 


ture (Goldberg & Bacon 

2005 Goldberg & Leonard 2007 

Okura & Futamase|2009 

Schneider & Er|2008 ). In our mass 


reconstruction we adopt the flexion formalism as defined 


by Cain et al. (2011 hereafter CSB). Our reconstruction 
method assumes measurements of the reduced flexion fields 

= J^/(I — /^) and T 3 = Q /{1 — from lensed source 
images using the Analytic Image Model (AIM) approach, 
rather than calculating flexion estimates from image char¬ 
acteristics like moments or shapelet weights. 

Using the estimators for each of the lensing fields ob¬ 
tained from analyzing cluster data images, we constrain a 
model lensing potential and the mass distribution which cor¬ 
responds to that potential. Our model is defined as lens po¬ 
tential values ^jjk on a “mesh” of points Ok to be constrained 
by the data. The value of the derivatives of the model po¬ 
tential at an arbitrary point 6 j are calculated using a gen¬ 
eralized finite differencing method using appropriate linear 
combination of nearby potential values on the mesh to give 
the desired derivative. The coefficients of this linear trans¬ 
formation are determined using a singular-value decomposi¬ 


tion (SVD) as in Bradac et al. (2009) to yield a system of 
equations relating the lensing potential values in the neigh¬ 
borhood of an arbitrary point in the field to the derivatives 
of the lensing potential: 


didTip{6j) = y^^ajk{n,m)'tpk, 


( 1 ) 
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where n and m are the derivative orders for each dimension, 
k indexes the mesh points and ajk{n,m) is the SVD matrix 
coefficient for the weight of the potential value contribution 
of the mesh point to the n, m derivative of the potential 
at Oj. As the value for each leasing field is a linear combi¬ 
nation of potential field derivatives, this scheme calculates 
each lensing field at any point in the field as a linear com¬ 
bination of model potential values from neighboring mesh 
points to Oj. To constrain the model at the location of a 
leased image, the model value for the lensing field(s) which 
have estimators derived from the image can be calculated 
and compared to the model. 

We then determine the values of which best fit the 
data by minimizing a model penalty function (x^) which 
includes contributions from weak lensing shear, strong lens¬ 
ing, and flexion to penalize models which do not produce 
the measured lensing distortions. The form of the strong 
and weak lensing penalty functions are identical to those 
of Bradac et al. (2009). The value of the flexion penalty 


function is determined by the reduced flexion fields and 
T 3 (where i indexes the Nf galaxies with measured flex¬ 
ion fields), and the model values for the flexion and conver¬ 
gence fields and calculated for the location of that 

galaxy: 
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( 2 ) 

Here, Np is the number of galaxy images with flexion mea¬ 
surements. As described above, the model flexion and con¬ 
vergence values are determined from a linear combination 
of mesh potential values. The uncertainties crjr,^ and ag^i 
include both intrinsic flexion scatter and measurement er¬ 
ror. This is the same basic form for the penalty function as 


described in Er et al. (2010), but we explicitly formulate it 
in terms of the reduced flexion estimators Ti and T 3 from 
CSB. 


As in Bradac et al. (2009), we include in our full penalty 
function regularization terms, which prevent the highly- 
nonlinear nature of the penalty function’s dependence on 
the mesh point potential values from driving the model to 
assume values wildly discrepant from the initial conditions 
during the reconstruction. This regularization is particularly 
important for mesh points in regions that are weakly con¬ 
strained by the data and therefore more susceptible to nu¬ 
merical effects in the reconstruction. We choose to regular¬ 
ize only on the values of convergence and shear, meaning 
that we penalize model potential values at each mesh point 
which require large pixel-to-pixel variations in the lensing 
fields fields relative to a reference value. Note that this does 
not mean that the model cannot vary from its initial condi¬ 
tions! Instead, we are requiring that the deviation from the 
initial conditions be driven by the data, rather than numer¬ 
ical effects. We iterate the reconstruction process and reset 
the reference value to the output of the previous iteration. 
This causes the model to step smoothly away from the ini¬ 
tial guess in a controlled and converging way. Without reg¬ 
ularization, the highly nonlinear dependence of the reduced 
lensing fields prevents the reconstruction process from find¬ 
ing a solution which minimizes the penalty function. 

Our regularization penalty function for convergence is 
defined over as a sum over the Nmesh mesh points located 


at Ok as 

^mesh 2 

xleg,K= XI 0 ^'^^ - K^^\0k)j , (3) 

k 

with similar terms for each component of the shear field. 
The value of the regularization weights and 

are determined empirically. We test a range of values for 
each reconstruction to not only find an appropriate value 
which prevents spurious nonlinear numerical artifacts while 
still allowing the reconstruction to deviate from initial con¬ 
ditions as indicated by data, but also to ensure that any 
reconstruction is stable with respect to variations the regu¬ 
larization coefficients. We also select values which place the 
total regularization penalty per mesh point (as a proxy for 
the regularization degrees of freedom) to be approximately 
equal that of the penalty per degree of freedom for each 
data portion of the full penalty function (i.e., strong, weak, 
or flexion lensing). 


3 TESTING THE RECONSTRUCTION 

METHOD 

With flexion included into the mass reconstruction, we test 
our method using simulated measurement catalogs. We gen¬ 
erate realistic lensing estimator values, including realistic 
measurement errors, from an underlying mass model which 
includes substructure. We then reconstruct the mass distri¬ 
bution using all the available data, and compare that recon¬ 
struction to a modified reconstruction which incorporates 
the same data for strong and weak lensing shear, but omits 
the flexion data. This allows us to make direct comparisons 
between reconstructions and evaluate the contribution that 
flexion makes to the fidelity of the final mass reconstructions. 
In particular, we are interested in finding substructure in the 
mass distribution which would be unobserved without the 
inclusion of flexion data. 

This latter point is important for how we construct our 
input simulated lensing measurements. It has already been 
well established that strong lensing is a sensitive probe of 
substructures. The goal of the mass reconstructions using 
our simulated measurement catalogs is to test whether re¬ 
alistic substructures in cluster lenses can exist which are 
detectable using flexion but not without. As we will show, 
this is indeed the case. 


3.1 Input Lensing Mass Distribution 


Our input mass distribution is created to mimic a galaxy- 
cluster sized halo with a small number of substructures, each 
simulated using a nonsingular isothermal ellipse model. We 
set the scale of the large central halo to be comparable with 
a galaxy cluster velocity dispersion of approximately 1500 
km/s. This is a large cluster, similar to those in the CLASH 
(Postman et al. | 2012 | l and Hubble Frontier Fields initiativ^ 
(PI Lotz: Program ID 13495) samples. We also include two 
substructures in the lens within a radius of 1.5 arcminutes 
of the cluster center. These subhaloes are also nonsingu¬ 
lar isothermal ellipses with varying masses from group scale 


^ http: / / www.stsci.edu/hst / campaigns / frontier-fields 
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down to large-galaxy scale (see Appendix A for tabulated 
subhalo masses). We construct six different data sets corre¬ 
sponding to substructured cluster haloes. 

The positioning of the halo mass substructures was not 
done in a fine-tuned way, i.e., we did not select the location 
or masses of the substructures in any way to enhance the 
effect of the substructures on the resultant lensing morphol¬ 
ogy. Instead, we use a conservative approach and restrict the 
possible locations such that the strong lensing morphology 
does not initially indicate that there is significant substruc¬ 
ture in the lens. The intent of this work is to characterize 
the differential effect that flexion makes on the resulting 
mass distribution. Thus, we restrict the substructure halo 
masses and positions to be outside the approximate radius 
of the critical curve so that the additional substructure is 
not easily detected using the strong lensing data alone. Sub¬ 
structures which are too massive or near the critical curves 
create significant deviations in the location, morphology, and 
magnification of strong lensing images. This can be an ex¬ 
tremely sensitive probe of substructure (e.g. Vegetti et al. 


2012 ), however our goal is to make a total census of the 


substructure in galaxy clusters. 


3.2 Simulated Lensing Estimators 


We simulate lensing measurement catalogs to test the mass 
reconstruction method and the effect of including flexion. 
Using the main halo and smaller substructures we define 
as above, we analytically calculate the appropriate lensing 
field estimators from the total mass profile using the full 
lens equation and its derivatives. We assume a single, mod¬ 
erate redshift of ziens = 0.3 for the lens halo (including its 
substructures) in all of the simulated datasets. Independent 
random noise is added to each simulated measurement with 
scatter chosen to match observational constraints from sin¬ 
gle orbit Hubble Space Telescope Advanced Camera for Sur¬ 
veys {HST/A.CS) cluster observations. For the simulations 
presented here, we employ an input lensing mass distribu¬ 
tion consisting of three non-singular isothermal ellipse (NIE) 
halo models for computational efficiency in constructing the 
simulated measurements and a simple comparison model. 

To determine our simulated measurement values based 
on the lensing mass distribution, the source-lens geometry 
must be determined along the line of sight in addition to 
the plane of the sky. This is because the lensing fields are 
all dependent on the redshift of the lensing mass as well as 
the redshift of the source object. We randomly select the 
redshift of each lensed source from a gamma distribution 

^ ^ exp(-V 2 o) (4) 


with zo = 2/3. This is the same distribution as used in 
previous work to emulate the true distribution of sources 


(Brainerd et al. 1996 Schrabback et al^|2007 Bradac et al. 


2004). We draw a new random redshift for each flexion and 


weak lensing source image, and each strong lensing image 
system has its own redshift as well. Redshift errors on the 
weak lensing/flexion sources are added with cr^ = 0.05(l+z). 
The strong lensing systems do not have redshift errors, to 
reflect the spectroscopic redshifts which are available for 
many multiply imaged galaxies. We describe the measure¬ 


ment generation procedure for each type of catalog (flexion, 
weak lensing, and strong lensing) below. 


To generate our flexion measurement catalog we select 
a random distribution of source positions for the flexion and 
leak lensing shear measurements. The positions are selected 
uniformly across a square field 3ffi on a side to a density 
of 80 galaxies/arcmin^. For the flexion measurements we 
calculate the reduced flexion fields and T 3 as defined in 
CSB at the position of each of the sources. We add Gaussian 
noise to each component of the flexion measurements with 
a standard deviation of 0.03 arcsec”^. These noise levels are 


motivated by previous studies of intrinsic flexion (Goldberg 
fc Leonard|2007[ GSB). 


We then remove from the flexion dataset any sources 
where the magnitude of the measured flexion exceeds 1.0 
arcsec“^. This ceiling on the measured flexion reflects how 
the quadratic expansion of the lens equation which defines 
the flexion fields is only valid when the flexion correction 
is small relative to the shear distortion. In a real dataset, 
images with flexion values in these regimes would either be 
identifiable as strongly-lensed arcs and thus in the strong 
lensing dataset instead of the flexion catalog, or would be 
removed during image fitting as in GSB. There the authors 
found that for flexion values this large the AIM image fit¬ 
ting did not produce accurate results, but the properties of 
the AIM covariance matrices allowed for efficient removal of 
these catastrophic outlier images. Our limit on the flexion 
values is a uniform approximation of the limit that would oc¬ 
cur in real data: the directly observable distortion in a lensed 
image due to flexion is quantified by the unitless product of 
the image scale (such as the half-light radius) and the re¬ 
duced flexion field, rather than the more physically relevant 
reduced flexion field alone. Because flexion is a quadratic 
effect, measuring flexion requires that the lensed galaxy im¬ 
age be larger than is required for ellipticity - in practice, 
this means 2-3 times the point-spread function size along 
the long axis of the image. In order to fully include varying 
image size and its effect on flexion estimation and the result¬ 
ing mass reconstruction, a better understanding of intrinsic 
flexion and its correlation to galaxy size and magnitude must 
be achieved. However, the simpler cut on the value of the 
flexion field is a sufficient approximation for our purpose. 


The weak lensing shear and strong lensing datasets are 
constructed from the reduced shear and deflections fields. 
The source locations for the shear are randomly distributed 
through the simulated field and locations with reduced 
shears nearing unity are removed (for the same reason as 
large flexion values are removed). The selection of the source 
position for the strong lensing systems is not completely ran¬ 
dom, however. For each system we initially draw a random 
position from a small field in the source plane centered on 
the main halo core and solve the lens equation to find the 
image positions. We then perform a visual inspection and 
remove systems if they show obvious signs of substructure. 
With some source positions relative to the distribution of 
mass in the lens, multiple image systems will form around 
the smaller substructure. We reject these systems, as the 
utility of strong lensing for substructure detection is already 
well established. Instead, we select only image systems that 
are representative of a single, central halo and do not have 
obvious substructure included. 
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3.3 Initial Conditions and Regularization 


Our approach to reconstructing the mass distribution from 
the simulated measurements is the same approach we would 
take using real data. We determine initial conditions for 
our model fitting, set a mesh of points where the model 
potential will be defined, and reconstruct the lensing po¬ 
tential all using data-driven conditions. We perform several 
reconstructions with different values for some of the less- 
strongly-constrained input parameters (e.g., the regulariza¬ 
tion weights). We evaluate the output reconstructions, and 
in particular we evaluate the behavior of the penalty func¬ 
tion values as the reconstruction converges, to determine the 
best values for these parameters. We do not see a strong de¬ 
pendence on the particular values - small changes do not 
significantly affect the output reconstruction. 

The first step in our reconstruction is to assign initial 
conditions for the reconstruction which are motivated by 
the data. We use an analytic model to set the initial model 
lensing potential values which, in turn, determine the initial 
value of the penalty function and its gradient. As was shown 
previously ( Bradac et al.pOOd Bradac et al.|20Q5 ), the final 
reconstruction using this method is not strongly dependent 
the particular initial model used provided that the regular¬ 
ization is chosen appropriately. We employ a simple, para¬ 
metric, axisymmetric model as our initial guess. We take a 
single, nonsingular isothermal sphere as our initial model, 
locating the center of the NIS at the average position of the 
images in the strong lensing data set. The mean of the dis¬ 
tances from these strong lensing images from their average 
position defines an approximate Einstein radius for the over¬ 
all lensing system. Using the approximated Einstein radius 
estimate we construct our initial model from a NIS pro¬ 
file. These simple, reasonable approximations provide data- 
driven initial conditions. 


Another important input parameter to the mass recon¬ 
struction is the weight of each component of the regulariza¬ 
tion in the penalty function etc., from Eq.[^, which de¬ 
fines the regularization weight relative to the strong lensing, 
weak lensing shear, and flexion contributions. As noted in 
^ too little weight given to regularization (small 77 values) 
yields reconstructions with unrealistically large amounts of 
small-scale spatial variation as well as aphysical lensing po¬ 
tential values, such as negative convergence. These numeri¬ 
cal artifacts of the nonlinear lens inversion prevent the opti¬ 
mization algorithm from finding the appropriate minimum 
as it iterates. Excessive weight to the regularization (large 77 ) 
causes the reconstruction to simply return a model which is 
minimally different from the initial mass model, preventing 
the data from informing the reconstructed mass distribution. 
Because it is not possible to determine the most appropriate 
77 values a priori, and because the most appropriate specific 
value for the regularization weight will depend not only on 
the details of the data available but also on the particular 
mesh geometry chosen, we determine 77 empirically by re¬ 
constructing with a range of values. By comparing the 
values for each part of the penalty function, we are able 
to separate successful reconstructions from those which are 
over- or under-regularized. 

By evaluating the values output by these reconstruc¬ 
tions we determine the approximate best combination of val¬ 
ues for aNis and regularization weight to use for the given 


dataset. To compare different reconstructions, we consider 
both the total x^ value, as an indicator of the overall fit¬ 
ness of the reconstruction, and the individual contributions 
of the different datasets. Reconstructions which yield values 
of x^/d.o.f.?^ 1 for each of the three datasets are preferred, 
even in cases where that increases both the overall and reg¬ 
ularization values. Eurthermore, we prefer models which 
reproduce the strong lensing data more precisely above those 
which minimize the weak lensing or flexion penalty func¬ 
tions, provided that there are no obvious signs of significant 
pixel to pixel variations. This is because strong lensing data 
have very small measurement errors, in fact smaller than the 
error we use in the strong lensing penalty function which we 
inflate in order to more easily locate the penalty function 
minimum. Thus if the strong lensing image systems are ro¬ 
bustly identified in a given cluster, we can reasonably require 
that the solution have a strong lensing penalty value of ef¬ 
fectively zero. The reconstruction which meets these criteria 
while minimizing the overall value is the reconstruction 
which we select as our solution. 

In the presence of a significantly substructured mass dis¬ 
tribution, the reconstruction will not converge with a single 
halo initial model, and will plateau to x^ values of ~5-10 (or 
more) per degree of freedom after several iterations. This is 
an indication of an initial guess which is far from the solu¬ 
tion, and the highly non-linear nature of the reconstruction 
process requires and update to the initial conditions based 
off this information. We update the initial conditions itera¬ 
tively, occasionally adding additional haloes at the location 
of and sizes of substructures detected in the non-converged 
reconstruction in order to arrive at a good fit. See 0 for a 
discussion of our substructure detection method. 

3.3.1 Reconstruction Model Mesh 

We distribute mesh points over the “observed field”, which 
in this case is the full area around the mean of the strong 
lensing image locations in a square region 3.^5 on a side (mim¬ 
icking the HST/ACS field of view). These mesh points are 
where the model lensing potential and its derived lensing 
fields are defined. The “base” mesh is defined by a rectan¬ 
gular, regular grid of points evenly defined across our field. 
We refine this grid in a circular region centered on the field 
center-point. The radius of this circular region is 1.5 times 
the mean distance of the strong lensing images from the av¬ 
erage strong lensing image location. Within this radius, we 
distribute mesh points in a rectilinear grid at twice the base 
density. 

We further refine in regions within a 3^^ radius of a 
strong lensing image position with four times the base den¬ 
sity. This refinement is necessary to accurately constrain the 
source positions for the strong lensing systems. Numerical 
error in the calculation of the gradient of the model lensing 
potential due to a paucity of mesh points near strong lens¬ 
ing images can dominate the strong lensing contribution to 
the penalty function, which will preclude convergence and 
accurate reconstruction. 

This nested mesh structure places increasing numbers 
of mesh points where the mass distribution is expected to 
be varying more rapidly, as well as where we have tighter 
constraints from the lensing data, and therefore a higher 
density of mesh points for the finite difference derivatives 
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used to calculate the model lensing field values enhances 
the accuracy of the reconstruction output. 

In the reconstructions presented here, we begin with 
a base mesh density of roughly 70 mesh points/arcmin^, 
corresponding to about one galaxy image per mesh point, 
with the density of mesh points higher in the refinement 
regions. The optimal configuration will have increased mesh 
point density wherever the potential is changing rapidly, and 
will not have any changes in refinement level of more than 
a single step. 


4 RESULTS 

Using the methods described above, we reconstruct the mass 
distribution for each of the simulated cluster datasets, and 
bootstrap the simulation 1 catalogs and reconstruct each 
of those. Here we present our results, first in a detailed in¬ 
vestigation of simulation 1 and the statistics of the recon¬ 
structions after bootstrap-resampling of the lensing catalogs, 
then more broadly for the rest of the cluster data sets. In 
discussing these results, we will refer to locations within the 
reconstructed field by celestial directions (north being up 
and east being to the left), as we would were these actual 
cluster fields and not arbitrarily oriented simulations. 


4.1 Simulation 1 


Figureshows a comparison of the input mass distribution, 
the reconstruction without flexion, and the reconstruction 
with flexion. In both the without-flexion reconstruction and 
in the with-flexion case, the strong lensing data is extremely 
well fit by the reconstruction, essentially perfectly repro¬ 
ducing the multiple image positions for each of the sources. 
The final weak lensing penalty function values in each re¬ 
construction are very similar as well. Both reconstructions 
fit the strong and weak lensing data equally well within sta¬ 
tistical uncertainty (we discuss the statistics of these recon¬ 
structions below). 

However, the reconstruction with flexion includes addi¬ 
tional substructure undetected in the flexion-less reconstruc¬ 
tion. One of the two subhaloes (to the northeast) is clearly 
detected by eye, and the second is apparent as an exten¬ 
sion to the southwest. To quantify substructure detections. 


we employ Source Extractor (SExtractor, Bertin & Arnouts 
1996) to detect “objects” in our mass maps. This is similar 


in some ways to detecting objects in a crowded field, though 
with an unusually fine pixel scale. We require a large min¬ 
imum area (an area of ~ 200 arcsec^) above the detection 
threshold to identify an “object”. This prevents SExtractor 
from identifying the small, local variations in the conver¬ 
gence map due to the mesh point spacings as individual 
objects. 

Because cluster subhaloes are not as discrete as galax¬ 
ies in a typical field and instead are a contiguous distri¬ 
bution of matter with significant blending between haloes, 
and because estimating the “background” convergence level 
in the environment of a parent halo is problematic, we do 
not discount morphological evidence for substructure in our 
reconstructions. While the southern substructure is not in¬ 
dividually detected by SExtractor, it instead appears as an 
elongation of the source in that direction. Combining both 



Figure 1. With flexion measurements we detect substructure in 
the mass distribution of Simulation 1. From top to bottom, the 
input convergence (surface mass density) map, the convergence 
map reconstructed without flexion, and the map reconstructed 
including flexion data. The innermost contour levels are the same 
in all three cases beginning at a convergence of 1.5 and decrease 
by 0.375 with each step, for a source at infinite distance. 
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the SExtractor detections and the visible morphology we 
conclude a strong detection of the northern subhalo and a 
likely detection of the second, though it blends partially into 
the main halo. 


4.2 Bootstrap Error Estimation 

We estimate the uncertainty in our mass reconstructions via 
a bootstrap resampling of the weak lensing shear and flex¬ 
ion catalogs. We randomly select measurements from the 
two catalogs uniformly, allowing entries to be selected mul¬ 
tiple times, to form a new pair of resampled catalogs. From 
these catalogs we perform a mass reconstruction using the 
same initial conditions model and reconstruction parame¬ 
ters as used in the second step reconstruction of the original 
dataset. We repeat this process 100 times, each with a dif¬ 
ferent realization of the resampled catalog, to generate a 
distribution of mass reconstructions which we then compare 
to the input convergence map. We hold the strong lensing 
sources constant in the bootstrapping process. This is ac¬ 
ceptable given that the identification of multiple images in 
real data allow us to more securely identify the real systems 
that the simulated strong lensing catalogs represent and the 
likelihood of mistakenly including interloper galaxies which 
are not strongly lensed is relatively low. 

The bootstrap reconstructions produce a large data- 
cube of information, with a distribution of model values for 
each of the lensing fields: -0, ai, k, 7 ^, and Qi. Ten in 
all, eleven if the magnification fi is included. To character¬ 
ize these distributions in a comprehensible way, we produce 
an average convergence map and the RMS variation of the 
convergence in the 100 reconstructions, which is displayed in 
Figure These maps show that the substructure detected 
with flexion is significant relative to the scale of statistical 
variations in the catalog data. 

We also map both the RMS variation of the boot¬ 
strapped convergence maps relative to the primary recon¬ 
struction, and also the difference between the primary and 
average convergence maps, presented in FigureThese give 
us a metric of bias in the primary reconstruction, i.e., how 
much of an outlier is the primary reconstruction. Were we 
to see a large amount of bias in these maps, meaning a large 
difference between the RMS variation about the average and 
the RMS variation about the primary, or were we to see a 
large difference between the primary and the average conver¬ 
gence maps, this would suggest that the primary reconstruc¬ 
tion were a statistical outlier rather than a robust minimum 
solution to reconstructing the mass distribution. We do not 
see any evidence that this is the case. 

We see the expected small number of outliers among 
the reconstructions, but overall the bootstrapped catalog re¬ 
constructions should not vary drastically from the primary 
reconstruction if they are valid estimators of the error in the 
reconstruction. A few reconstructions (^5) have hugely dis¬ 
crepant penalty function values, universally due to a poor 
fit to the the strong lensing data and thus an incorrectly 
located critical curve. These outliers are easily removed and 
discarded from our statistical analyses, though leaving them 
in does not significantly affect the results as the lensing field 
maps, such as the convergence, are not catastrophically de¬ 
viant from the primary solution even with these outlier so¬ 


lutions since they depend less on the exact position of the 
critical curve. 

The reconstruction method is computationally inten¬ 
sive, requiring the inversion of large sparse matrices used 
for the finite-differences calculation of the gradient of the 
penalty function with respect to the lensing potential values 
on the mesh, so we perform this bootstrap analysis on only 
one of our simulations (Simulation 1 ), both with and with¬ 
out flexion measurements included, to estimate the amount 
of variation to expect due to the statistical uncertainty in 
our lensing field estimator measurements. As with real data, 
the variation depends not only on the accuracy and density 
of the weak lensing and flexion estimators across the field, 
but also on the amount of specific constraint provided by the 
particular set of strong lensing images available. The num¬ 
ber and, importantly, the location of the multiple images 
have an effect on the spatially varying amount of freedom 
that the mass model can change within the statistical errors 
of the measurements. A more full statistical analysis of the 
cluster lensing mass reconstruction made including flexion 
allowing a broader variation of possible halo and subhalo 
configurations is a large undertaking, and beyond the scope 
of this work. However we can use our bootstraps to begin 
to estimate the errors in our mass maps, and compare the 
difference between the input mass distribution and the re¬ 
constructed mass map to these error estimates. 

Using the bootstrap-resampled catalog reconstructions, 
we can understand the statistical variation in our mass maps 
due to the catalog sampling. One immediate feature to note 
in the non-flexion RMS map (Figure]^ is that the variation 
is small and is relatively uniform across the field. This means 
that despite the catalog resampling, a consistent mass dis¬ 
tribution is preferred by the data. The input subhaloes are 
larger than any of the statistical variations by quite a large 
margin, the non-detection of the subhaloes is not due to 
noise or a statistical fluctuation. Instead, it is due to incom¬ 
plete information in the lensing inputs, namely the lack of 
flexion information measuring the mass and shear gradients. 

The RMS map for the flexion reconstruction is less 
uniform than the non-flexion RMS map, which can be un¬ 
derstood considering how flexion measurements inform the 
mass reconstruction. Because the flexion held is sourced very 
locally by the gradient of the mass distribution, and since the 
hexion measurement is not dominated by statistical noise in 
the same way that shear measurements are, excluding cer¬ 
tain individual hexion measurements from the catalog can 
have a more signihcant effect on the resulting mass map 
than happens with weak lensing shear measurements that, 
by their nature, must be correlated together to extract the 
lensing signal. Thus we expect that bootstrapping the hex¬ 
ion measurement catalogs in addition to the shear catalogs 
will increase the RMS values near where the hexion mea¬ 
surements are most informative and constraining. This in¬ 
cludes the locations of substructures, as can be seen in the 
RMS maps in Figures [2|3| but also where the hexion held is 
strong but reliable measurements are relatively sparse (near 
the main halo center) as well as where hexion measurements 
are particularly near a strong lensing image. The RMS map 
shows that neither substructure detection can be attributed 
to the statistics of the hexion information, with one cer¬ 
tain detection (north) and a second, less signihcant detec¬ 
tion (south). 
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Figure 2. Flexion-detected substructure is robust to catalog resampling. On the left, the average convergence map for the set of 
reconstructions using the bootstrapped catalogs (above) and the RMS variation about that average (below) for reconstructions of data 
from simulation 1 without flexion. On the right, the same two plots for simulation 1 reconstructions including flexion. Contours for the 
top two plots are as in Figure^ For the bottom two plots, the innermost contours begin at 0.2 and decrease by 0.025 with each step, 
also assuming a source at inflnite distance. 


As noted above in comparing the average reconstruc¬ 
tion to the primary reconstruction, both with and without 
flexion, we see no significant overall bias (Figure]^. For re¬ 
constructions with flexion, the average reconstruction under¬ 
estimates the mass in the substructures as is expected since 
the bootstrap reconstructions include catalogs where some 
of the flexion measurements nearest the substructures are 
not included which will reduce the lensing signal from that 
substructure. The features in these maps are what would be 
expected given the information available (and not available) 
to each reconstruction. 

These results, the primary mass reconstruction together 
with the properties of the reconstructions from bootstrap 
resampled catalogs, show that the inclusion of flexion mea¬ 
surements can detect substructures which are otherwise un¬ 
detectable in strong lensing and weak lensing shear data 
alone. Furthermore, these detections are statistically signif¬ 
icant and robust to the variations of bootstrapped catalogs. 


i.e., they do not depend solely on a few fortuitously-located 
images, though there is a more marked dependence of the 
resulting mass map on the location of the flexion measure¬ 
ments relative to any substructures than is seen for weak 
lensing shear measurements. 

4.3 Simulations 2 through 6 

Figures and each show comparisons between the re¬ 
constructed mass distributions (on the left) and the input 
mass distribution (right) the eight simulations. These, along 
with Table [AT] show the varying substructure sensitivity en¬ 
hancement that including flexion yields. Simulations 1 and 
2 each have two equal mass subhaloes, though the positions 
differ between the simulations. For simulations 3-6, we vary 
the masses of the subhaloes. 

A few trends emerge from these simulations. As would 
be expected, more massive substructures are reliably de- 
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Figure 3. The RMS variation of the bootstraps about the primary reconstruction (top) and the difference between the primary and 
the average bootstrapped reconstructions (bottom) for simulation 1. The left column is without ffexion, the right with ffexion. Green 
contours indicate positive values, with the innermost contours at 0.2 decreasing in steps of 0.025 to zero. Blue contours indicate negative 
values, with the innermost value at -0.2 and the same step size. 


tected, particularly if they are well separated from the core 
of the main halo. The one instance where a larger subhalo 
is not detected (simulation 7) is one where the subhalo is 
blended into the main halo. In many configurations, this 
type of subhalo (in terms of mass and location) would be 
detectable from strong lensing data alone, as we know from 
the set of halo+substructure configurations we generated 
but rejected for this study due to their obviously substruc¬ 
tured strong lensing observations. 

For smaller substructures, the detectability decreased 
with distance from the main halo. This is also expected. The 
value of the non-reduced flexion fields (i.e., T and 5) is dom¬ 
inated by the gradients in the convergence and shear field 
whose distortions depend on the proximity to the subhalo 
mass distribution. The total convergence (which factors into 
the reduced flexion) is dominated by the main halo mass. 
At larger radii from the main halo center, the total conver¬ 
gence falls away from unity, the flexion decreases quickly 
away from the subhalo, and both effects cause the measur¬ 


able reduced flexion signal to decrease. The least massive 
subhalo we detect (in simulation 4), though we only detect 
it marginally, was located quite close to the main halo center 
where the total convergence is much closer to unity. More ra¬ 
dially distant subhalos were either undetected or more mas¬ 
sive. As we discuss in detail in §4.4| the inclusion of flexion 
does more than detect substructure and improve the recon¬ 
structed ellipticity. Flexion also improves the overall fidelity 
of the mass reconstruction independent of the detectability 
of individual structures. 


4.4 Aperture Masses 

In addition to detestability of substructure in terms of iden¬ 
tifying individual haloes, we can also quantify the fidelity 
of our mass reconstructions in terms of how how well the 
mass maps compare to the input mass distribution by com¬ 
paring aperture masses. Figure plots these results. We 
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Figure 4. A comparison of input mass distributions (left) and the reconstructions including flexion (right) for simulations 1 (top) and 
2 (bottom). Contours are as in FigureIn these plots, and in Figures [5|6| the reconstruction held displayed is centered on the strong 
lensing image centroid, and therefore are slightly offset from the true halo center. 


compare aperture masses rather than point-by-point masses 
since we expect relationship between the input and recon¬ 
structed mass distribution to not only include noise propa¬ 
gated from the lensing field estimators, but also an effective 
convolution kernel that will vary somewhat across the field 
based on the density of lensing field measurements and the 
type of measurements in that region (i.e., strong lensing, 
shear, or flexion). 

There is a moderate tightening of the scatter (in dex) as 
the aperture radius increases, which is consistent with the 
idea of the reconstruction method imposing and effective 
convolution kernel just as larger apertures in photometry 
will omit fewer photons dispersed by a non-trivial point- 
spread function. Though the measurements for a single halo 
location with multiple aperture radii are certainly corre¬ 
lated, there is similarly good agreement between the input 
mass distribution and the reconstructed masses in each of 
these apertures, confirming that we are in fact reconstruct¬ 
ing the substructures on these small scales, as the fraction of 


the mass enclosed due to the substructure increases as the 
aperture shrinks. 

The measurements scatter about the unity line by 13% 
for reconstructions with flexion, which together with the 
RMS scatter between the different bootstrapped catalog re¬ 
alizations which have an RMS variation of ~ 10% suggest an 
intrinsic scatter of approximately 8%. Without flexion, on 
the other hand, the reconstructed aperture masses are sys¬ 
tematically biased low at the location of the substructures, 
as we expect given that the substructures are not detected 
without flexion, by 25-40%. This underestimate is also seen 
in aperture measurements of the main halo masses, show¬ 
ing that the improvement from adding flexion is not only a 
substructure phenomenon. 


4.5 Simulation Summary 

We see from these simulations the feasibility for detecting 
and identifying substructures depends both on their mass 
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Figure 5. As in Figure^ for simulations 3 (top) and 4 (bottom). 


and position relative to the cluster halo. We do not have 
many detections of substructure at small radii because we se¬ 
lect “un-substructured” strong lensing image conhgurations 
and are therefore artificially enhancing the likelihood that 
substructure near the Einstein radius of the cluster will re¬ 
main undetected. More centrally located subhaloes are also 
less likely to be detected because of the reduced number 
of flexion measurements — again, these substructures are 
those that would typically be visible in the strong lensing 
data but do not appear in our data because we select against 
the configurations that make them apparent. However, in¬ 
dependent of subhalo identification, the masses we measure 
in apertures as small as 10^^ are accurate to an error (scatter 
plus statistical) of 13% across the reconstructed field. This 
is much smaller than the observed systematic offset from not 
including flexion in the reconstruction. 

Flexion data included into the mass reconstructions fill 
an important information gap in the full lensing data set. 
With a mass limit that depends both on the angular size 
of the subhalo Einstein radius (which is a proxy for the 
subhalo mass) and its distance from the main cluster halo 


center, flexion can resolve cluster substructures which are 
otherwise unobservable from the strong and weak lensing 
data alone. Though the specihc approach to appropriately 
identify an individual subhalo in these mass reconstructions 
has room for optimization, the aperture mass results show 
that the accuracy of the reconstruction is high, and that re¬ 
constructions without flexion systematically underestimate 
these aperture masses. This is true even away from the sub¬ 
halo locations, where flexion information better constrains 
the overall main halo shape as well. 


The three-dimensional halo mass of the substructures 
at the margins of our detection threshold in this simulation 
sample have masses of 2-3x10^^ M© within an aperture of 
10", though there is a dependence on the radial position 
of the subhalo in its detectability. This demonstrates the 
efficacy of flexion as a probe of small-scale galaxy cluster 
substructure. 
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Figure 6. As in Figure]^ for simulations 5 (top) and 6 (bottom). 


5 DISCUSSION 

We have shown that including flexion into the mass recon¬ 
struction enhances the sensitivity of the mass reconstruc¬ 
tion to substructure a galaxy-cluster-scale lens. Significantly 
massive subhaloes which otherwise are undetected in data 
sets including only multiple image systems and ellipticity 
measurements are made detectable by including flexion data 
measured to a precision achievable using current techniques 
(e.g., CSB). The additional information from flexion requires 
no additional observational investment - single-orbit HST 
observations are sufficient for including flexion into the lens- 
ing analysis and the mass reconstruction for a typical cluster. 

This result indicates that the addition of flexion into a 
lensing mass reconstruction better constrains the formation 
and structure of galaxy clusters and the subhalo mass func¬ 
tion. The work here shows the distinct possibility that dark 
subhaloes in galaxy clusters of significant mass, such as those 
included in the simulated measurements we present here, are 
not detected in lensing mass maps produced without flex¬ 
ion. Constraints on the subhalo mass function depend on the 
number density and mass of the detected subhaloes within a 


cluster, as well as a thorough understanding of the detection 
limits for subhaloes of a given mass. By not including flexion, 
the number density of detected substructures will only be a 
lower limit on the true substructure density. With such large 
structures (up to roughly 10% the mass of the main halo) it 
is unlikely, though not impossible, that a cluster would have 
that much mass in substructure ( Giocoli et al.|20'T0 ). Given 
that substructures this large were not detectable without 
flexion, to constrain smaller amounts of substructure with 
confidence requires that flexion data be included. 


The possibility of having significant substructure in a 
galaxy cluster field which is otherwise undetected by strong 
or weak lensing analyses also has important implications on 
the inferences drawn from the lensing mass distributions. 
For example, recent work using galaxy clusters as “cosmic 
telescopes” to select and study high-redshift (^ ^ 7) galaxies 


(e.g. Hall et al.|2011 Postman et al.|20T2 Zheng et al.|2012 


[Zitrin et al. 2012) and constraining the properties galaxy 
population formed at or shortly after the epoch of reioniza¬ 
tion requires an accurate magnification map of each cluster 
field. The overall cluster halo amplifies the effect of sub- 
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Figure 7. The input aperture mass is accurately recovered by 
the flexion reconstructions. Color indicate different aperture radii, 
with blue, green, and red for 10", 15", and 20", respectively. The 
line indicates unity, not a fit to the data. Points marked with an 
X are measured from reconstructions without flexion, and sys¬ 
tematically underestimate the input mass. The RMS variation of 
bootstrapped catalog reconstructions indicate a typical aperture 
mass error of 10%. 


structures on the resulting magnification map, as much of 
the cluster is typically very near the leasing critical density. 
By including flexion into the reconstruction and more ac¬ 
curately constraining the leasing potential, the systematic 
error in the inferred intrinsic properties of the high-redshift 
galaxy population from unresolved substructures is reduced. 

Using Simulation 1 as an example, as compared to an 
identical main halo without substructure, a HST/ACS field 
containing the substructured halo lens would probe only 
~ 88% of the solid angle that would be inferred to have been 
probed if the single halo were the only one reconstructed. 
Though the magnitude of this systematic error in each case 
is dependent on the exact position and mass of the substruc¬ 
tures, it is not unreasonable to estimate that for any clusters 
analyzed without flexion, there is a 10-15% systematic un¬ 
certainty in the solid angle, and therefore the differential 
volume, probed at any given source redshift. 

Furthermore, the intrinsic luminosity determined for 
any detected high-redshift objects is dependent on the local 
absolute-value of the magnification determined by the lens 
reconstruction. Again referring to the comparison between 
our simulated lens with substructure and the halo without 
substructure, the average absolute magnification across a 
HST/ACS field for an infinitely distant source is ~ 50% 
higher for the substructured lens, which would correspond 
to an average shift of 0.45 magnitude. Locally the variations 
between the two lens models is much larger, and in some re¬ 
gions will magnify while in other demagnify a source if the 
substructures are detected. 

The distribution of local variations in absolute magni¬ 
fication has a long tail corresponding to deviations in the 
critical curves which would most likely only affect strongly 
lensed/multiply imaged sources in the z > 7 population. In 


the rest of the field, where a large number of the observed 
high-redshift galaxy sources are likely to be observed, the 
ratio of the substructured lens magnification to the single¬ 
halo lens is well-correlated to the substructured lens mag¬ 
nification (p ^ 0.75). This means that in regions where the 
magnification is more enhanced by the substructure, i.e., 
where magnification maps without flexion will more likely 
be erroneous, non-detection of substructure and the system¬ 
atic error in the magnification map produced by that non¬ 
detection will create a larger error in the inferred intrinsic 
luminosity of the object. Typical magnification ratios in re¬ 
gions near the substructures are \psub/fJ^singiel ~ 5, which 
would create a shift of ~ 1.75 magnitude in the intrinsic 
source magnitude inferred without substructure. That said, 
there is a significant spread to be expected in these magnifi¬ 
cation ratios based on the image position relative to the crit¬ 
ical curve. The Einstein radius of the substructure imbedded 
in the main halo in our simulations vary from about 
meaning that individual detections can be even more than 
the typical factor and, for lensed image positions near the 
centers larger substructures the inferred magnification could 
in fact be lower than expected from a non-substructures 
lens model. This type of systematic error would significantly 
change the resulting luminosity function, as well as the in¬ 
ferred properties of the lensed sources. 


Lensing can also determine cluster mass distributions to 
be compared with other mass estimators. In particular, the 
Sunyaev-Zel’dovich decrement and the X-ray surface bright¬ 
ness of a cluster, as functions of position, also constrain 
the mass distribution of the cluster under the assumption 
of hydrostatic equilibrium for the intracluster medium. N- 
body simulations show that non-thermal pressure support, 
from cosmic rays and bulk gas motion, plays a significant 
role in the pressure budget of the cluster, accounting for 


up to half of the total pressure supporting the cluster (Nel¬ 


son et al. 2014). While these effects can be accounted for 


statistically for large cluster samples, doing so on the level 
of substructures within the cluster requires additional in¬ 
formation. Cluster mass profiles, with the substructure ap¬ 
propriately constrained using flexion, provide an essential 
test for the hydrostatic assumption, allowing the fractional 
non-thermal pressure support to be directly quantified for 
individual clusters. 


Flexion is an important addition to the variety of lens¬ 
ing information available from detailed imaging data. The 
results presented here strongly motivate the application of 
this reconstruction technique to simulated data sets with 
more realistic mass distributions than the simple halo pro¬ 
files we employ, and additionally applying it to real data. 
And while there remains important work to be done to make 
flexion as robust as the more mature lensing techniques 
(strong lensing and weak lensing shear), e.g. improving our 
measurement of intrinsic shape scatter in flexion measure¬ 
ments, better understanding the effects of different source 
selection on the flexion noise properties, etc., we have shown 
that this work is well worth the undertaking for the resulting 
enhanced accuracy of reconstructed mass and magnification 
maps that come from including flexion. 
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APPENDIX A: EXTENDED APERTURE MASS 
RESULTS 

Table |A1| shows input and measured masses for each sub¬ 
halo, along with the distance from each subhalo to the clus¬ 
ter halo center. We include a three-dimensional mass for the 
subhalo alone, as well as the mass within the aperture for 
the full input mass distribution, and the reconstructed mass 
distribution. Note that the aperture masses include contri¬ 
bution from both the subhalo and the main halo. Values are 
tabulated at different aperture radii - 10'', 15", and 20". We 
also include aperture mass measurements for the locations 
of the subhaloes in Simulation 1 as reconstructed without 
ff exion. 
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masses. * denotes reconstructions without flexion) 


Cluster substructure with gravitational flexion 
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